#include "IO1D.h"

#include <sstream>

#include "../Input/PicParams.h"
#include "../Parallel/ParallelMPI1D.h"
#include "../Field/Field1D.h"
#include "../PhysicalField/PhysicalField.h"
#include "../Species/Species.h"


using namespace std;

IO1D::IO1D(PicParams* params_, ParallelMPI* pmpi_) :
    IO(params_, pmpi_)
{


}

IO1D::~IO1D()
{
}

//create hdf5 data hierarchical structure: datespace, dateset and so on
void IO1D::createFieldsGroup()
{
    fieldsGroup.dims_global[2] = params->dims_space_global[0] + 1;
    fieldsGroup.dims_global[1] = 1;
    fieldsGroup.dims_global[0] = 1;

    fieldsGroup.ndims_[0] = fieldsGroup.dims_global[0];
    fieldsGroup.ndims_[1] = fieldsGroup.dims_global[1];
    fieldsGroup.ndims_[2] = fieldsGroup.dims_global[2];


    fieldsGroup.offset[0] = 0;
    fieldsGroup.offset[1] = 0;
    fieldsGroup.offset[2] = 0;


    fieldsGroup.stride[0] = 1;
    fieldsGroup.stride[1] = 1;
    fieldsGroup.stride[2] = 1;


    fieldsGroup.block[0] = 1;
    fieldsGroup.block[1] = 1;
    fieldsGroup.block[2] = 1;


    //for attribute
    fieldsGroup.aDims = 3;

    //============ Create fieldsGroup ============================================

    addField(fields->rho_global);
    addField(fields->phi_global);
    addField(fields->Ex_global);
    addField(fields->rho_global_avg);
    addField(fields->phi_global_avg);
    addField(fields->Ex_global_avg);
    for(int i = 0; i < fields->rho_s.size(); i++)
    {
        addField(fields->rho_s_global[i]);
        addField(fields->rho_s_global_avg[i]);
        addField(fields->Vx_s_global_avg[i]);
        addField(fields->Vy_s_global_avg[i]);
        addField(fields->Vz_s_global_avg[i]);
        addField(fields->Vp_s_global_avg[i]);
        addField(fields->T_s_global[i]);
        addField(fields->T_s_global_avg[i]);
        addField(fields->q_s_global_avg[i]);

        if(params->is_temperature_anisotropic)
        {
            addField(fields->Tx_s_global[i]);
            addField(fields->Tx_s_global_avg[i]);

            addField(fields->Ty_s_global[i]);
            addField(fields->Ty_s_global_avg[i]);

            addField(fields->Tz_s_global[i]);
            addField(fields->Tz_s_global_avg[i]);
        }

    }
    fieldsGroup.dataset_id.resize( fieldsGroup.dataset_stringName.size() );

}


//Create particles h5 file pattern
void IO1D::createDiagsGroup()
{
 
}



//write potential, rho and so on into hdf5 file every some timesteps
void IO1D::write( PicParams* params_, ParallelMPI* pmpi_, PhysicalField* fields, vector<Species*>& vecSpecies, Diagnostic* diag, int itime)
{
    const char* h5_name;
    int n_dims_data = 3;
    int iDiag;

    hid_t       group_id;
    hid_t       dataspace_id;
    hid_t       memspace_id;
    hid_t       dataset_id;
    herr_t      status;

    hsize_t     dims_global[3];
    hsize_t     count[3];             //size of subset in the file
    hsize_t     offset[3];            //subset offset in the file
    hsize_t     stride[3];
    hsize_t     block[3];

    //set stride and block as default
    stride[0] = 1;
    stride[1] = 1;
    stride[2] = 1;

    block[0] = 1;
    block[1] = 1;
    block[2] = 1;


    Diagnostic1D* diag1D = static_cast<Diagnostic1D*>(diag);


    if(itime % params->timestep_io == 0 && pmpi->isMaster())
    {
        string step_output_string;

        step_output = itime / params->timestep_io;
        step_output_string = to_string(step_output);
        
        if(step_output_max >= 10 && step_output_max <100)
        {
            if(step_output < 10)
            {
                step_output_string = "0" + step_output_string;
            }
        }
        else if(step_output_max >= 100 && step_output_max <1000)
        {
            if(step_output < 10)
            {
                step_output_string = "00" + step_output_string;
            }
            else if(step_output<100)
            {
                step_output_string = "0" + step_output_string;
            }
        }
        else if(step_output_max >= 1000)
        {
            log_warning<<"step_output_max is too large, please change the code in IO1D.cpp";
        }

        data_file_name = "data/data" + step_output_string + ".h5";
        data_file_id = H5Fcreate( data_file_name.c_str(), H5F_ACC_TRUNC, H5P_DEFAULT, H5P_DEFAULT);

        //============= write attributes, n_dim: number of dimension ======================
        hid_t attrs_dataspace_id, attrs_id;
        int n_dim = 1;
        hsize_t attrs_dims[1];
        attrs_dims[0] = 1;
        attrs_dataspace_id = H5Screate_simple(1, attrs_dims, NULL);
        attrs_id           = H5Acreate2(data_file_id, "n_dim", H5T_STD_I32BE, attrs_dataspace_id, H5P_DEFAULT, H5P_DEFAULT);
        H5Awrite(attrs_id, H5T_NATIVE_INT, &n_dim);
        H5Sclose(attrs_dataspace_id);
        H5Aclose(attrs_id);

        //=============write fields============================================
        fieldsGroup.group_id = H5Gcreate(data_file_id, "/Fields", H5P_DEFAULT, H5P_DEFAULT,H5P_DEFAULT);
        for(int i = 0; i < fieldsGroup.dataset_stringName.size(); i++)
        {
            fieldsGroup.dataspace_id = H5Screate_simple(n_dims_data, fieldsGroup.dims_global, NULL);
            h5_name = fieldsGroup.dataset_stringName[i].c_str();
            fieldsGroup.dataset_id[i] = H5Dcreate2(fieldsGroup.group_id, h5_name, H5T_NATIVE_DOUBLE, fieldsGroup.dataspace_id,
                                                    H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            fieldsGroup.status = H5Dwrite(fieldsGroup.dataset_id[i], H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, fieldsGroup.dataset_data[i]);
            fieldsGroup.status = H5Sclose(fieldsGroup.dataspace_id);
            fieldsGroup.status = H5Dclose(fieldsGroup.dataset_id[i]);
        }
        fieldsGroup.status = H5Gclose(fieldsGroup.group_id);


        //============== write Diagnostic: particle_flux, heat_flux and angle_distribution ============
        //============== particle_number, particle_energy, radiative_energy_collision      ============

        group_id = H5Gcreate(data_file_id, "/Diagnostic", H5P_DEFAULT, H5P_DEFAULT,H5P_DEFAULT);

        //total_electric_field_energy
        dims_global[0] = 1;
        dims_global[1] = 1;
        dims_global[2] = 1;

        dataspace_id = H5Screate_simple(n_dims_data, dims_global, NULL);
 
        dataset_id   = H5Dcreate2(group_id, "total_electric_field_energy", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(diag1D->total_electric_field_energy));
        status       = H5Dclose(dataset_id);
        status       = H5Sclose(dataspace_id);

        //particle_number, total_particle_energy, particle_flux_left, particle_flux_right, heat_flux_left, heat_flux_right
        dims_global[0] = 1;
        dims_global[1] = 1;
        dims_global[2] = diag1D->n_species;

        dataspace_id = H5Screate_simple(n_dims_data, dims_global, NULL);
 
        dataset_id   = H5Dcreate2(group_id, "particle_number", H5T_NATIVE_INT, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_INT, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(diag1D->particle_number[0]));
        status       = H5Dclose(dataset_id);

        dataset_id   = H5Dcreate2(group_id, "total_particle_energy", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(diag1D->total_particle_energy[0]));
        status       = H5Dclose(dataset_id);

        dataset_id   = H5Dcreate2(group_id, "particle_flux_left", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(diag1D->particle_flux_left[0]));
        status       = H5Dclose(dataset_id);

        dataset_id   = H5Dcreate2(group_id, "particle_flux_right", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(diag1D->particle_flux_right[0]));
        status       = H5Dclose(dataset_id);

        dataset_id   = H5Dcreate2(group_id, "heat_flux_left", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(diag1D->heat_flux_left[0]));
        status       = H5Dclose(dataset_id);

        dataset_id   = H5Dcreate2(group_id, "heat_flux_right", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(diag1D->heat_flux_right[0]));
        status       = H5Dclose(dataset_id);

        status       = H5Sclose(dataspace_id);    

        //angle_distribution
        dims_global[0] = 1;
        dims_global[1] = diag1D->n_species;
        dims_global[2] = 90;

        dataspace_id = H5Screate_simple(n_dims_data, dims_global, NULL);

        dataset_id = H5Dcreate2(group_id, "angle_distribution_left", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        for(int ispec = 0; ispec < diag1D->n_species; ispec++)
        {
            offset[0] = 0;
            offset[1] = ispec;
            offset[2] = 0;
            count[0]  = 1;
            count[1]  = 1;
            count[2]  = 90;

            memspace_id = H5Screate_simple(n_dims_data, count, NULL);
            status      = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, offset, stride, count, block);
            status      = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, memspace_id, dataspace_id, H5P_DEFAULT, &(diag1D->angle_distribution_left[ispec][0])) ;
            status      = H5Sclose(memspace_id);
        }
        status = H5Dclose(dataset_id);

        dataset_id = H5Dcreate2(group_id, "angle_distribution_right", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        for(int ispec = 0; ispec < diag1D->n_species; ispec++)
        {
            offset[0] = 0;
            offset[1] = ispec;
            offset[2] = 0;
            count[0]  = 1;
            count[1]  = 1;
            count[2]  = 90;

            memspace_id = H5Screate_simple(n_dims_data, count, NULL);
            status      = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, offset, stride, count, block);
            status      = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, memspace_id, dataspace_id, H5P_DEFAULT, &(diag1D->angle_distribution_right[ispec][0]));
            status      = H5Sclose(memspace_id);
        }
        status = H5Dclose(dataset_id);
        status = H5Sclose(dataspace_id);

        //energy_distribution
        dims_global[0] = 1;
        dims_global[1] = diag1D->n_species;
        dims_global[2] = diag1D->n_energy;

        dataspace_id = H5Screate_simple(n_dims_data, dims_global, NULL);

        dataset_id = H5Dcreate2(group_id, "energy_distribution_left", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        for(int ispec = 0; ispec < diag1D->n_species; ispec++)
        {
            offset[0] = 0;
            offset[1] = ispec;
            offset[2] = 0;
            count[0]  = 1;
            count[1]  = 1;
            count[2]  = diag1D->n_energy;

            memspace_id = H5Screate_simple(n_dims_data, count, NULL);
            status      = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, offset, stride, count, block);
            status      = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, memspace_id, dataspace_id, H5P_DEFAULT, &(diag1D->energy_distribution_left[ispec][0]));
            status      = H5Sclose(memspace_id);
        }
        status = H5Dclose(dataset_id);

        dataset_id = H5Dcreate2(group_id, "energy_distribution_right", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        for(int ispec = 0; ispec < diag1D->n_species; ispec++)
        {
            offset[0] = 0;
            offset[1] = ispec;
            offset[2] = 0;
            count[0]  = 1;
            count[1]  = 1;
            count[2]  = diag1D->n_energy;

            memspace_id = H5Screate_simple(n_dims_data, count, NULL);
            status      = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, offset, stride, count, block);
            status      = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, memspace_id, dataspace_id, H5P_DEFAULT, &(diag1D->energy_distribution_right[ispec][0]));
            status      = H5Sclose(memspace_id);
        }
        status = H5Dclose(dataset_id);
        status = H5Sclose(dataspace_id);

        //psi_rate
        dims_global[0] = 1;
        dims_global[1] = 1;
        dims_global[2] = diag1D->n_psi;

        dataspace_id = H5Screate_simple(n_dims_data, dims_global, NULL);
 
        dataset_id   = H5Dcreate2(group_id, "psi_rate_left", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(diag1D->psi_rate_left[0]));
        status       = H5Dclose(dataset_id);

        dataset_id   = H5Dcreate2(group_id, "psi_rate_right", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(diag1D->psi_rate_right[0]));
        status       = H5Dclose(dataset_id);
        status       = H5Sclose(dataspace_id);

        //radiative_energy_collision
        dims_global[0] = 1;
        dims_global[1] = diag1D->n_collision;
        dims_global[2] = diag1D->params->dims_space_global[0];

        dataspace_id = H5Screate_simple(n_dims_data, dims_global, NULL);

        dataset_id = H5Dcreate2(group_id, "radiative_energy_collision", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        for(int i_collision = 0; i_collision < diag1D->n_collision; i_collision++)
        {
            offset[0] = 0;
            offset[1] = i_collision;
            offset[2] = 0;
            count[0]  = 1;
            count[1]  = 1;
            count[2]  = diag1D->params->dims_space_global[0];

            memspace_id = H5Screate_simple(n_dims_data, count, NULL);
            status      = H5Sselect_hyperslab(dataspace_id, H5S_SELECT_SET, offset, stride, count, block);
            status      = H5Dwrite(dataset_id, H5T_NATIVE_DOUBLE, memspace_id, dataspace_id, H5P_DEFAULT, &(diag1D->radiative_energy_collision[i_collision][0]));
            status      = H5Sclose(memspace_id);
        }
        status = H5Dclose(dataset_id);
        status = H5Sclose(dataspace_id);

        //sigma_left, sigma_right
        dims_global[0] = 1;
        dims_global[1] = 1;
        dims_global[2] = diag1D->n_species;

        dataspace_id = H5Screate_simple(n_dims_data, dims_global, NULL);

        dataset_id   = H5Dcreate2(group_id, "sigma_left", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(diag1D->sigma_left[0]));
        status       = H5Dclose(dataset_id);

        dataset_id   = H5Dcreate2(group_id, "sigma_right", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(diag1D->sigma_right[0]));
        status       = H5Dclose(dataset_id);

        status       = H5Sclose(dataspace_id); 


        //vdf, edf
        if(diag1D->params->diagnostic_param.size() == 1)
        {
            dims_global[0] = 1;
            dims_global[1] = diag1D->diagnostic_param.vedf_n_species;
            dims_global[2] = diag1D->diagnostic_param.vedf_n_velocity;
            dataspace_id = H5Screate_simple(n_dims_data, dims_global, NULL);

            dataset_id   = H5Dcreate2(group_id, "vdf_velocity_array", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, diag1D->vdf_velocity_array.data_);
            status       = H5Dclose(dataset_id);

            dataset_id   = H5Dcreate2(group_id, "edf_energy_array", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, diag1D->edf_energy_array.data_);
            status       = H5Dclose(dataset_id);


            dims_global[0] = diag1D->diagnostic_param.vedf_n_position;
            dims_global[1] = diag1D->diagnostic_param.vedf_species_number_list.size();
            dims_global[2] = diag1D->diagnostic_param.vedf_n_velocity;
            dataspace_id = H5Screate_simple(n_dims_data, dims_global, NULL);

            dataset_id   = H5Dcreate2(group_id, "vdf_parallel_to_B", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, diag1D->vdf_parallel_to_B.data_);
            status       = H5Dclose(dataset_id);

            dataset_id   = H5Dcreate2(group_id, "vdf_perpendicular_to_B", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, diag1D->vdf_perpendicular_to_B.data_);
            status       = H5Dclose(dataset_id);

            dataset_id   = H5Dcreate2(group_id, "vdf_z", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, diag1D->vdf_z.data_);
            status       = H5Dclose(dataset_id);

            dataset_id   = H5Dcreate2(group_id, "edf", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
            status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, diag1D->edf.data_);
            status       = H5Dclose(dataset_id);
            //close dignostic group
            status = H5Gclose(group_id);
        }

        //============== write ParticleSource load 1d: load_dn, load_number_rem_tot, load_source_density, load_source_density_pre ============
        group_id = H5Gcreate(data_file_id, "/ParticleSource", H5P_DEFAULT, H5P_DEFAULT,H5P_DEFAULT);

        int n_particle_source = vecParticleSource->size();
        vector<double> val_temp;
        val_temp.resize(n_particle_source);
        dims_global[0] = 1;
        dims_global[1] = 1;
        dims_global[2] = n_particle_source;

        dataspace_id = H5Screate_simple(n_dims_data, dims_global, NULL);

        //load_dn
        for(int i = 0; i < n_particle_source; i++)
        {
            if((*vecParticleSource)[i]->ps_param.ps_type == "load")
            {
                auto ps1d = static_cast<ParticleSource1D_Load*>((*vecParticleSource)[i]);
                val_temp[i] = ps1d->load_dn;
                //cout<<"load_dn: "<<ps1d->load_dn<<endl;
            }
            else
            {
                val_temp[i] = 0.0;
            }
        }
        dataset_id   = H5Dcreate2(group_id, "load_dn", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(val_temp[0]));
        status       = H5Dclose(dataset_id);

        //load_number_rem_tot
        for(int i = 0; i < n_particle_source; i++)
        {
            if((*vecParticleSource)[i]->ps_param.ps_type == "load")
            {
                auto ps1d = static_cast<ParticleSource1D_Load*>((*vecParticleSource)[i]);
                val_temp[i] = ps1d->load_number_rem_tot;
            }
            else
            {
                val_temp[i] = 0.0;
            }
        }
        dataset_id   = H5Dcreate2(group_id, "load_number_rem_tot", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(val_temp[0]));
        status       = H5Dclose(dataset_id);

        //load_source_density
        for(int i = 0; i < n_particle_source; i++)
        {
            if((*vecParticleSource)[i]->ps_param.ps_type == "load")
            {
                auto ps1d = static_cast<ParticleSource1D_Load*>((*vecParticleSource)[i]);
                val_temp[i] = ps1d->load_source_density;
            }
            else
            {
                val_temp[i] = 0.0;
            }
        }
        dataset_id   = H5Dcreate2(group_id, "load_source_density", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(val_temp[0]));
        status       = H5Dclose(dataset_id);

        //load_source_density_pre
        for(int i = 0; i < n_particle_source; i++)
        {
            if((*vecParticleSource)[i]->ps_param.ps_type == "load")
            {
                auto ps1d = static_cast<ParticleSource1D_Load*>((*vecParticleSource)[i]);
                val_temp[i] = ps1d->load_source_density_pre;
            }
            else
            {
                val_temp[i] = 0.0;
            }
        }
        dataset_id   = H5Dcreate2(group_id, "load_source_density_pre", H5T_NATIVE_DOUBLE, dataspace_id, H5P_DEFAULT, H5P_DEFAULT, H5P_DEFAULT);
        status       = H5Dwrite (dataset_id, H5T_NATIVE_DOUBLE, H5S_ALL, H5S_ALL, H5P_DEFAULT, &(val_temp[0]));
        status       = H5Dclose(dataset_id);


        status       = H5Sclose(dataspace_id);

        //close particle source group
        status = H5Gclose(group_id);

        //=========================== close file ======================================================
        //close hdf5 file
        status = H5Fclose(data_file_id);




    }


}
